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Abstract 

Granular Solid Hydrodynamics (GSH) is a broad-ranged continual mechanical description of 
granular media capable of accounting for static stress distributions, yield phenomena, propagation 
and damping of elastic waves, the critical state, shear band, and fast dense flow. An important input 
of GSH is an expression for the elastic energy needed to deform the grains. The original expression, 
though useful and simple, has some draw-backs. Therefore, a slightly more complicated expression 
is proposed here that eliminates three of them: (1) The maximal angle at which an inclined layer 
of grains remains stable is increased from 26° to the more realistic value of 30°. (2) Depending 
on direction and polarization, transverse elastic waves are known to propagate at slightly different 
velocities. The old expression neglects these differences, the new one successfully reproduces them. 
(3) Most importantly, the old expression contains only the Drucker-Prager yield surface. The new 
one contains in addition those named after Coulomb, Lade-Duncan and Matsuoka-Nakai - realizing 
each, and interpolating between them, by shifting a single scalar parameter. 

PACS numbers: 46.05.+b ; 43.25. +y; 62.20.D; 45.70.Cc; 
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I. INTRODUCTION 



GSH (brief for: Granular Solid Hydrodynamics) is a continual mechanical theory 



constructed to account for a broad range of granu 



distribution 



5[, yield 



-ML incremental stress-strain relation 
of elastic waves 7]], elasto-plastic motion jg], the critical state 



ar phenomena, including static stress 



propagation and damping 
9j, shear band and fast dense 
flow [lOj. An important input of GSH is an expression for the elastic energy needed to 



deform the grains 



111 ] . The energy density w(uij) is a function of the elastic strain Uij, 



which we define as the long-scaled, coarse-grained measure of how and how much the grains 
are deformed. Therefore, both the elastic energy w and the stress (which comes about 
only because the grains are deformed) are necessarily functions of u^. Moreover, we have 



o~ij = —dw/du 



i.i ■ 



(1) 



because is closely related to the force, and to the coordinate. So is given if w is. 

The relation between the elastic and total strain ey is only simple at small increments, 
where Suij ~ Se^ holds. More generally, it is given by the evolution equation of 0, 0, Is], 
and not by a function, because only state variable, not e^. As it turns out, the 

critical state is simply the stationary solution of Ui/s evolution equation [9]. 

The expression w(uij) is an input of GSH. It may either be deduced microscopically, 
or obtained iteratively in a trial and error process, by comparing the ramification of the 
proposed expression with experimental observations. As the first method is notoriously 
difficult, we choose trial and error. The original expression, 



w 



BVA (2A 2 /5 + u 2 Ji) , 



(2) 



"u"u- K 



is very simple, and a function of only two invariants, A = — Ukk and u 2 s 

5ij/3 is the traceless or deviatoric strain.) B and ^ are two (density dependent) 
elastic coefficients. The dependence on the third invariant v% = u*k u kj u ji * s neglected. 
Nevertheless, a number of granular features are contained in this expression. First of all, the 

n Q 

measured velocity of elastic waves [12| are well rendered [7j. Second, satisfactory agreement 
was achieved Q with the incremental stress-strain relation as observed and reported in [3] . 
In both cases, one is looking at small increment of the elastic strain, 8v,ij ~ fey, and the 



calculation employing Eq fl2]) is purely elastic (or hyperelastic): 

x d(Ji i x d%w x ro\ 

d<J ij = q du k£ = -~ n OEM- {o) 

dun duijduke 

Third, the Drucker-Prager yield surface, requiring a granular system at rest to have a ra- 
tio of shear stress a s = ^Jo*~a\- and pressure P = a#/3 smaller than a certain value, 
a s jP < constant, is, as explained next, an integral part of this expression. 

Generally speaking, in a space spanned by stress components, there is a surface that 
divides two regions in any granular media, one in which the grains necessarily move, another 
in which they may be at rest. This surface is usually referred to as the yield surface. Aiming 
to make its definition more precise, we take the yield surface to be the divide between one 
region in which elastic solutions are stable, and another in which they are not - clearly, the 
medium may be at rest for a given stress only if an appropriate elastic solution is stable. 
Since the elastic energy of a solution satisfying the equilibrium condition VjOj = is an 
extremum lj], the elastic energy is convex and minimal in the stable region, concave and 
maximal in the unstable one. In the latter case, an infinitesimal perturbation suffices to 
destroy the solution. The elastic energy of Eq (J2]) is convex only for 



u 



a /A < ^2B/A, implying a s /P A < ^2A/B. (4) 



In this paper, we propose a slight generalization of the energy, by including the third 
invariant u| = u^u^u*^ and its elastic coefficient %■> as 

W = B^a( 2 -A^ + U-^4). (5) 



.5 T £ A 
We shall in the rest of the paper examine its ramifications. 

There are many different yield surfaces in the literature. The first was proposed by 
Coulomb, who observed that slop es of sand piles never exceed a critical value and saw an 



analogy with the friction law 15], 



; =sm<p, (6) 

where ip is the internal friction angle, a material parameter, and a\ < 02 < 03 are eigenvalues 



of cry, ordered by their magnitude. Although the Cou 



estimating the stability of granular materials at rest 16j, a number of other models are also 



omb yield model is widely used for 



frequently used by engineers, depending on personal preference and experience, especially 



those by Drucker and Prager[l3|, Lade and Duncan 18l |. Matsuoka and Nakai 19|, given 
respectively as 



p 



3 + sin ip 



< J i (J 2< J 3 (1 — sin ip) cos 2 ip 
27P 3 ~ 



(3 — sin ipY 
(cti - o 2 f 



8 tan ip. 



(7) 
(8) 
(9) 



CTiCr 3 a 2 <73 <J\<J 2 

As we shall see below, all three yield models are also reproduced in satisfactory approxima- 
tion by the new energy expression. In addition, it also improves on other residual discrepan- 
cies between the elastic theory and measurements, such as in the incremental stress-strain 
data, or the speed of elastic waves. All this indicates that the new energy, now with three 
parameters: B,£,x, is capable of giving a more accurate description of granular elasticity. 

We note that the new cubic term in Eq (J5J) was first introduced by Humrickhouse in 14j, 
in an attempt to increase the maximum angle # max of inclination, at which a granular layer 
remains stable. Unfortunately, he only considered negative values of x- As these did not 
yield any improvement, he abandoned this term. As we shall see in section II, a positive x 
does yield a larger 6* max , of around 30°. In section III, we deduce the yield surface associated 
with Eq (JHJ); in section IV, incremental stress-strain relation and granular acoustics are 
studied. All support a positive x- Section V contains discussion and conclusions. 



II. THE MAXIMUM ANGLE OF INCLINATION 



We consider noncohensive granular materials, with constant mass density p. Denoting 

Al = U XXJ A 2 = Uyy, A3 = U ZZ , 



A4 U X y, A5 U X Z1 Ag Uy Z , 



the 6x6 Hessian matrix of the function w(A a ) is 

d 2 w 



H, 



(10) 



(11) 



8A Q dA p ' 

with a, {3 = 1, 2, 6 and the eigenvalues h\ < h 2 < ■ ■ ■ < h 6 . The associated yield surface, 
written as a function of the stress components, is given by 



h x = 0. 
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(12) 



As the eigenvalues frequently lack analytic expressions, we consider instead its determinant, 
det H a/ 3 = 0. This is of course only a necessary condition, as we need to ensure that the 
vanishing eigenvalue is the smallest one, hi, while the other ones are larger. This is done 
numerically 

Next, consider an infinite granular layer in gravity, inclined by an angle 9, see FigJTJ The 



elastic strain is taken to assume the form 
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Uij = A 



r 

t -1 



(13) 



with r = u xz /A, where r 2 = u 2 J (2 A 2 ) - 1/3. Inserting the strain Eq ([13]) into (IA2IIA6D . 



-0.31 




FIG. 1. Variation of the maximum angle of inclination with the parameters £ and x calculated 
with the potential Eq ([5]) and the strain Eq (|13p. 
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we obtain 



v., 



o 7 



2 + x 3/2 

-A T . 



9£ + 5 X + 15 A 3 /2 x + 2 3/2 2 
9£ 2£ 



(14) 
(15) 



In addition, because of the force balance V,<7jj = pg^ the two stresses are also related by 



c TZ = &zz tan 9. 



(16) 



It is worth remarking that if we insert the strain Eq ( TT3|) into the Hessian matrix H a p, we 
find for x > that two eigenvalues are 4 — x (l ± 3 y/l + 4 r 2 ) , implying that a stable layer 
exists only for x < l, because at least one eigenvalue is negative for x > 1. As the other 
eigenvalues cannot be expressed analytically we now consider det (H a p) = 0, finding 



with 



det(H a p) = A{2 + x) A 3 D 1 D 2 , 



D 1= (-4 + 2x + 2x 2 + 9xV), 

D 2 = 27x 2 (x + 2) r 4 + 3 (-30* + 12 X 2 + 5 X 3 + 36£ X 2 + 12) r 2 
+2(x-l) (4 X + X 2 -6)+36e (X 2 -1)- 



(17) 



(18) 
(19) 



The real roots of D\ = are t±= ±^4 — 2x 2 — 2x/3x, that of D 2 = 



arc 



with 



■3xV 2(2 + x)' 
h = -5 x 3 - 12 (3f + 1) x 2 + 30 x - 12, 



X 



72^x 5 + 12 (I08e 2 - 5) x 4 - 24 (72£ + 11) X 3 



(20) 

(21) 
(22) 



+36 (48£ + 25) x 2 - 720x + 144. 



(Because ki < for x G (0, 1), positive square root for \fk~2~ ie taken here.) Requiring the 
vanishing eigenvalue to be the smallest one, we find yield to occur at r = r c , implying a 
maximum angle of inclination by employing Eqs f Tl4|T5"]) . 

18 (2 + x)r c 



9 V , 



=p arctan 



L9(x + 2)t c 2 + 18£ + 30 + 10x. 



(23) 
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Inserting r c of Eq (|20|) into it, and remembering that the sign only indicates a left or right 
inclination, we obtain the absolute value for the maximal angle as 



^max — arctan 



6 y / 2xv / 2Txv / £:i + Vh 
h + Vh + 36 £x 2 + 60x 2 + 20 x 3 



(24) 



It is nice to have an analytic expression for # max , which reduces to the results of H] for 



X -> i.e.: r c -> ^/f - 1/3 and # c -> arctan [y9f -3/ (3f + 2)] . 

In figUl the expression of Eq ([24"]) for 6* max is plotted against £ for various < x < 1- 
The biggest value # max = 30° is achieved at £ ~ 0.71 and x — 0.31. More specifically, as x 
increases from to 0.31, the peak of 6* max versus £ increases too, but decreases after that. 

A final remark: the elastic strain of Eq ( fl~3|) is a result of assuming that displacement 
vector Ui varies only with the layer depth z, and has nonvanishing components along the x 
and z directions, U y = 0. Then yield, at 9 = 8 m8iX , occurs simultaneously in the whole layer. 
This is an idealization. In reality, the displacement vector should also have a nonvanishing 
y component. Then yield will probably start at the layer top. This case will be studied 
elsewhere. 



III. THE YIELD SURFACE 

In this section, we consider uniform stress cx^ and density p, and employ the coordinate 
system of the principle directions of the stress. The elastic strain is also uniform and in 
its principle system, ie, a^- = and Uij = for i ^ j, see Appendix. 

Instead the principle strains u xx , u yy , u zz , we shall, for simplicity, use A, S > and L, 
where 



U 



vv 



u 7 



A 

T 

A 

~3 
A 



1 — S sin f L 



7T 



3/J 



1 + SsmL) 



1 -5 sin IL + ^ 



(25) 
(26) 
(27) 



(In soil mechanics L is usually called as Lode angle, of which value range is taken as [0, 2%)). 
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And we have the following expressions for the principle stresses, 

(Xr-r _ n: S 



3V3 + — (l2cosL-4V3 sin L + V3S^j (28) 



A 3 / 2 4£ 



(6\/3 cos2L — 18 sin 2L + V3S sin 3L 



°yy 



xs 2 
+- — 

24£ 

A s/2 -3V3+^S(S + 8smL) (29) 

_ \/3x 5 2 / 12 CQS 2L-S sin 3L) , 
24£ v ; 

= 3^3 - — (l2 cos £ + 4^3 sin L- \/3<s) (30) 



f 18 sin 2L + 6^ cos 2L + v^S sin 3L 
24£ V 

The pressure and the deviatoric stresses q, Q are simpler, with A factored out, 

P = a kk /3 = ^ A 3 / 2 U + \S 2 + sin 3Lj , 



(31) 



q = a zz - a xx = — A 3/2 5 ( x 5 sin 2L - 4 cos L) , (32) 
\/3 

Q/3 = (P - <Jyy) = ^—A 3/2 S (xS cos 2L - 2 sin L) . (33) 



A. The Cylindrically Symmetric Case 

For a xx = tjyy and the sample is cylindrically symmetric, implying a Lode angle 

L = 71 /Q or 77r/6, see Eqs. (125|26j) . Cylindrical symmetry is usually assumed in analyzing 
the so-called "triaxial test," widely used in soil mechanics. Inserting these Lode angles into 
Eqs (I31|32p . we obtain 

q_ 18 S f(x«S-4) forL = vr/6, 

P xS s + 65 2 + 72^ X I ( X<S + 4) for L = 7vr/6. 

Note L = 7r/6 is the case of triaxial extension with q < 0, while L = 7tt/6 is the case of 
triaxial compression with q > 0. At yield, 5 = S y i e u, Eq ( IMj) implies q ~ P, and in the 
stress space spanned by P, g yields two straight lines. 

To obtain the value for S yie id, we calculate the determinant det(H a/3 ), obtaining 

27A 3 



det(tf a/3 ) = ^P 2 P 2 P 3 , (35) 
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and for L = 77r/6: 



where for L = 7r/6: 

Di = X S+2, (36) 

D 2 = X S-A, (37) 

D 3 = X 2 S 4 - 8 X S 3 - 2AS 2 - UAZxS + 288£; (38) 

D l = - X S+2, (39) 
D 2 = - X S-A, (40) 
D 3 = X 2 S A + 8 X S 3 - 2AS 2 + 144£xS + 288£. (41) 

If x — 0, the equation det(H a p) = reduces to D 3 = 48S 2 — 576£ = for both Lode angles, 
implying S y i e id = \/12£ for extension and compression. Inserting this into Eq fl34l . we have 
q/P = ±a/3/£, the two straight yield lines are then symmetric about the P-axis. 

For x 7^ 0, one needs to resort to numerical calculation, which shows that S y i e id ^ is 
given by the smallest positive root of D 3 = 0. Yet because Eqs (1381) and (jHJ) are different, 
Syieid is different for extension and compression, and so are the two slopes, see figf5]-a. This 
unsymmetric behavior of yield is well documented and familiar in soil mechanics. Within 
the present framework, it is a measure of how much the third invariant influences granular 
elasticity. Figf2]-b shows variation of the slope of yield line with the parameter £, at varying 
^ % < 1. In agreement with the observations, the absolute value of the slopes is always 
greater for compression than that for extension, see figj2]-c. 

B. The General Case 

Next we consider the general case, with a xx ^ a yy ^ o zz . Inserting the expressions 
(I25|26ll27p into the Hessian matrix, and calculating its determinant, we have 

27A 3 

det(^) = j^rDiD2, (42) 

where 

D 1 = x 3 S 3 sin3L-6x 2 S 2 + 32, (43) 
D 2 = X 3 S 5 sin 3L - 6 X 2 S A - 40x^ 3 sin 3L (44) 
-48 (3£x 2 + l)<S 2 + 576£. 
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FIG. 2. Yield behavior as determined by the energy Eq ([5|), for the case of cylindrical symmetry, 
(a) Straight yield line for £ = 1 and x = 0-2. (b) Variation of the slope of yield line with £ for 
X = 0,0.2,0.4,0.6,0.8,1 (increasing indicated by arrows), (c) The difference between the slopes for 
compression and extension. 



It can be demonstrated numerically that for any L G [0, 2n), yield S = S y i e id is given by the 
smallest positive root of Di 2 = 0, which is also computed numerically. Inserting the yield 
strain (L, S yie id) into (I28|I29|I30I) the yield surface in the stress space can be plotted. Since 
the compressional strain A is an overall factor, the surface is conelike, as illustrated in figJH 
We present the yield surface (as is customary in soil mechanics) with a closed curve given by 
its intersection with the so called 7r-plane, defined by P =const., see figj3) For the conelike 
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FIG. 3. Yield surface, 7r-plane and their intersection loci in principle-stress space. The 7r-plane is 
defined by constant a xx + a yy + a zz . 

surface it is convenient to introduce the rectangular coordinates in the 7r-plane defined by 

7Ti = [a zz - o yy ) /y/2P, (45) 
vr 2 = (2a xx - a yy - a zz ) / V6P. (46) 

Using Eqs (I28|I29|I30I) . they become 
6V6S (Sx 

sin2L — 4cosL) 

1X1 = x(sin3L)S 3 + 6<S 2 + 72£ ' ( ' 

[3 6S 2 x cos 2L — 24S sin L + 24a/3<S cos L — 6a/35 2 x sin 2L 
712 ~\ 2 X (sm3L)S 3 + 6S 2 + 72£ ' [ } 

Inserting (L, S y i e id) into them the yield curve in the 7r-plane can be plotted, see figJH On the 

other hand, the yield surfaces of the Coulomb, Drucker-Prager, Lade-Duncan and Matsuoka- 

Nakai models, also conelike, and their loci in the 7r-plane, may be plotted directly - using 

Eqs (I6]r7|8|9|) and Eqs ( I45|46|) - and compared. 

The Drucker-Prager yield model can be reproduced exactly by the energy Eq ()5]) with 

X = 0. It is a circle in the 7r-plane, see figJU-a. As x increases from 0, the yield curve 

of Eq (|5]) is distorted from a circle to a triangle-like curve, in fairly good agreement with 

the Lade-Duncan and Matsuoka-Nakai models, see figJl]-b,c. Although the Coulomb model 

is a hexagon, with three more corners, see the dotted line of figJU-b, the difference seems 

insignificant. 
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FIG. 4. Yield loci in the 7r-plane calculated with the energy Eq ([5]), in comparison with various 
yield models, from (a) to (d): Drucker-Prager, Coulomb, Lade-Duncan, Matsuoka-Nakai. 

We note the large difference between the maximum angle of inclination 8 max and the 
friction angle ip. Employing the parameters of Fig.6 below: £ = 1, x = 0.25, we obtain 
#max = 29.8°, and ip = 65°, 54°, 60° for the Coulomb, Lade-Duncan, and Matsuoka-Nakai 
model, respectively. 



IV. THE COMPLIANCE TENSOR 



Since GSH is a unified description of granular behavior, the energy of Eq (jSJ) is expected 
to also account for the elastic stiffness and the speed of elastic waves, both experimentally 
accessible. The stiffness tensor Mij mn is of fourth order and defined as 

8(Tij = M ijmn 5u mn = ~(d 2 W I 'dUijdUmnjSUmn, (49) 
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FIG. 5. Comparison of compliance coefficients between those measured in [13f] and calculated via 
Eq ([MD. 



while the compliance tensor is its inverse, 



dUij X.ij mn d,0~ rnn . 



(50) 



The differential forms of Eqs P9|50l) are also referred to as incremental stress-strain relation 
in soil mechanics. Physically, the components of the stiffness tensor can be interpreted as 
response coefficients of a small stress change to a strain increment, or vice versa for those of 
the compliance tensor. 

Systematic measurements of the response coefficients as functions of stress and density 
was carried out by Kuwano and Jardine for the case of cylindrical symmetry using a triaxial 
apparatus 13J. A comparison of their data with the GSH calculation obtained for x — was 
given in j^la]. Here we consider whether the agreement may be further improved by a finite 
X- In figEJ we plotted the four compliance components, denoted as E xx = —l/X X xxx, E zz = 
—t/^zzzz, Gxz = — {/4\ xzxz and G xy = —l/^\ X yxy, as functions of P. The agreement is 
clearly improved for x — 0.1. Especially noteworthy is the fact that G xz and G xy degenerate 
for x — 0, but are split appropriately, in the correct order, for x = 0.1, again indicating that 
the sign of x is positive. 
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V. ELASTIC WAVES 



A further important elastic property is the speed q of elastic waves, which can be calcu- 
lated from the square root of the appropriate eigenvalues of the acoustic tensor, see [3, I20I ]. 



q k m k n d 2 w (51) 

13 pk 2 du im du nj 



where k m is the wave vector of propagation, with k 2 = k m k m . Inserting the energy of Eq (jSJ) 
into (I5ip . and eliminating the strain with the help of Eq (JT]), we can obtain different velocities 
as functions of the stress: q = c, (<7y). Such a calculation with x — has been given in 
resulting in satisfactory agreement, though the weak variation of the transverse wave velocity 
with varying polarization and direction of propagation was not reproduced. Here, we focus 
on this and the additional difference taking ^ / 0. Since the energy is a homogenous 
function of the strain, or equivalently the stress, the stress dependence of the velocities can 
be factorized as c$ ~ p 1 / 6 f [k m , £, x, a^JP) , where the factor / ~ c^/P 1 / 6 represents the 
influences of shear stresses on sound speeds (from which the values of the parameters £, x 
may be obtained with great accuracy). Fig. |6] shows the calculated q/P 1//6 , setting £ = 1 and 



X = 0.25, in comparison with the data reported in [12|. Cylindrical symmetry is assumed, 
and the wave is taken to propagate along axial (Fig. O-a) and radial directions (Fig. [B]-b). 
Then Cj/P 1//6 depends only on q/P, or equivalently on q/a xx = 3(g/P)(3 — q/P)^ 1 . Clearly, 
GSH reproduces the data fairly well, where especially the order of splitting of the transverse 
waves velocity, when either polarized along the axial or the radial direction is correct, once 
again indicating that x > 0. 

Finally, we note that irrespective of the energy expression, there can only be two different 
transverse velocities in a cylindrical geometry along the principle axes. Therefore, if hard- 
ened experimental evidences to the contrary arises, this would indeed be a reason to include 
an additional variable such as the intrinsic anisotropy. 

Assuming a general expression for the energy, dw = . . . dA + . . . du s + . . . dut, with 
dA = —5ijduij, du s = u^jduij/ug, du t = Ui m u m j/u 2 , we calculate the stress as 

&ij = A (A, u s , utjdij + A(A, u s , u t )uij + A(A, u s , u t )u im u mj , (52) 
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FIG. 6. Comparison of sound 
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speeds as calculated employing the energy expression of Eq ([5]) with 
12|. Pq is the atmosphere pressure. 



where fi depend on the energy w. Similarly, the stiffness tensor Mijki is generally given as 



M, 



ijkl — 9l8ijfikl + 92 {Uijdkl + UklSij) + g 3 (u im U m j5kl + Uk m U m i5ij) 



+ 94 (Sikdji + 5u5jk) + g^UijUki + g 6 {uijU km u ml + u k iu im u mj ) 

+ g7U im UmjUkmU m i + g s (ukiSij + UkjSu + uudjk + uijSik) (53) 

with gi = gi(A,u s ,u t ) again to be calculated from w. In the system of principle axes, 



U: 



U(i)6ij, the stiffness tensor reads 



ijkl 



9i + 92 + «(fe)) + 9hU{i)U{k) + 5-3 (tifi) + W( fc )) 
+#6 (w(fc)W(i) + u 2 (i)U(k)) + 97U 2 {i) u 2 {k) dijdki 

9a + 9s {u(i) + U(j)) {SikSji + SuSjk) 
g{i, k)5ij5ki + h(i,j) (8 ik 5ji + 5u8j k ) ■ 



(54) 
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The acoustic tensor Cu ~ Mijkirijnk (with rii = ki/\k\) therefore is 



Cu ~ (g(i,l) + h(i,l))mni 

+ (94 + 9s (u(i) + [u 2 + (it! - u 2 ) n\ + (u 3 - u 2 ) nl])) 5 U 

= H(i, VjUiUi + (#4 + 5f 8 (U( i} + u dev )) 6u, 

with H(i, I) = g(i, I) + h(i, I) and Ud ev (n) = u 2 + (ui — u 2 ) n\ + (u 3 — u 2 ) n\. 

For cylindrical symmetry, denoting z as the preferred direction, u\ = u 2 7^ w 3 , we have 
H(l, 1) = H(2,2) = H(l,2) and Ud ev — ui + (u 3 — ui)n\. We have, if the wave vector is 
along z (713 = 1, n\ — n 2 — 0), 



/ 



C z 



94 + 9s («i + W3) 

94 + 98 (ui + u 3 ) 

H 33 + g 4 + 2g 8 u 3 J 

if it is along x {n\ — 1, n 2 — n 3 — 0), 



V 



\ 



H n +94 + 2g 8Ul 

C x . = g± + 2g%u\ 

y g 4 + g 8 («3 + «i) y 

if it is along y (ra 2 = 1, ni = n 3 = 0), 



/ 



(7^ — 



\ 



# 4 + 2g 8 U! 

H u +g 4 + 2g 8Ul 
y g 4 + g 8 («3 + mi) y 

We note two different eigenvalues of the acoustic tensor for the transverse case, 



C 2 ,ll — Cz,22 — Cx,33 — C?/,33 7^ Cx,22 — Cj/,11, 



(55) 



(56) 



(57) 



(58) 



and realize that there can only be two different velocities. Note also that without g 8 that 
stems from the third invariant, only one transverse velocitiy exists. 



VI. DISCUSSION AND CONCLUSIONS 



Macroscopically speaking, granular media at rest are elastic, and one can account for all 
its static, mechanic behavior including yield, compliance coefficients, and sound propagation 
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by a single potential, the elastic energy. Finding an quantitatively appropriate potential is 
certainly interesting and useful, and the expression given in Eq ([5]) seems a good starting 
point. It is a generalization of the potential given in {4, 11], with a term depending on 
the cubic strain invariant. All empirical yield models and the other experiments considered 
above support this term. As mentioned in the introduction, yield models in soil mechanics 
are usually selected by personal experience, preference or convenience, and there seems no 
consensus which one is best. The uncertainty of the situation may be reduced with the 
help of the considered potential, as it unifies these models, and link them to other easily 
measurable elastic properties such as sound speeds, and compliance tensor. In this context, 
we would like to stress the importance of more accurate and systematic experiments, such 
as simultaneous measurement of sound velocity and yield. 



It is important to realize that yield as considered here happens at the highest possible 
stress at which the system may maintain an elastic, static solution. It is different from 
that obtained under stead shear: When approaching the critical state from an isotropic one, 
the approach is not monotonous if the starting density is high, and the limiting, so-called 
critical stress is lower than some of the stress values the system has undergone. At all these 
stress states, if the strain rate is stopped, the system will stop as well and retain its stress 
statically. Therefore, it is still below yield at all these stress states. Measurements carried 
out with steady shearing samples therefore do not reveal the yield surface. Interpreting the 
results as such will lead to discrepancies. 



Moreover, density also influences elastic properties of granular materials, which was con- 
sidered in [l|, but is neglected here for simplicity. In jl|, we have shown how to include the 
density dependence of the energy w to account for the so-called "cap" [2l|, or the sound 
velocity as measured by Hardin and Richart 22| . These features are easily transferred to the 
energy of Eq (jHJ). Summarizing, we believe that although the details of granular elasticity 
is complicated, and more accurate experiments are needed for further clarification, Eq (jSJ) 
provides a reasonable start point. 
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Appendix A: Stress-Strain Relation 

From the energy of Eq (jSj), we may calculate the stress as a function of the elastic strain 
via cry = —dw/duij, obtaining 

_ °K, y u xz u yz ^xy^zz) \X T ^ ) J-iU'xy ( \T\ 

&xy = 7t^- ' l A -V 

(A2) 



3X 


[U XZ Uy Z 


^XyUzzJ 


(X + 2) Au^ 










3x 


(u x yUy Z 


^xz^yy) 


(X + 2) Am m 










3x 


(U x y1lxz 


UyzUxx) 


(x + 2) A M?/2 



= 7S? ' (A3) 



and 



1 _ X /a"/ \ I 3X U yz U xy U xx + W L / A ,s 

q = o zz -a xx = — VA(ii„-!J 22 ) + - = , (A4) 

S/^ £ VA 

1 - X nr / x 3 x M L - u lz ~ u lx + u lv < a — \ 

q = a yy -a xx = -^-VA{u xx -Uyy)+-r— ^= -■ ( A 5) 

9e + 5 X + 15 A3/2 3( X + 2) 

= ~ A 1 H — v A (u xx + u yy ) (A6) 

+ («L + UxxUyy + u\ y + u 2 xy + u| 2 + uj z ) 

+ ( u xx + 4 ^*% + u 2 yy - 2u 2 xy + u 2 xz + u 2 yz ) 



+ 



2£\/A v "' a:a! ' ~"" xj: "" yv ' ""ff ""at/ 1 1 

2A 3 / 2 ^ ^ Wa; 2/ _ U xx u yy) ( u yy + ^ra) — U xz U yy ~ U yz U xx + ^UxyUxzUyz] 



Clearly, if the strain is diagonal, the stress is too. In the principle coordinate, the stress- 
strain relations becomes 

1 - X fir i. . \ i 3x ~u 2 xx + u\ 



<1 



u2 VA(u xx -u zz ) + ^ xx ^ z \ (AT) 



and 



9 = 772^ A («» - «w) +f ^ ( A8 ) 



9£ + 5 X + 15 a3/2 3(x + 2) ^ 

= A^ H — V A {Uxx + Uyy) (A9) 

+ y/~A£ { U XX + U XxUyy + 



+ ( M L + ^xxUyy + U 2 yy ) 

~^ 2A 3 / 2 ^ ^ _ UxxUyy ) 
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